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FOREWORD 


This  report  addresses  one  aspect  of  the  one-dimensional  surface  heat-flux  calculations 
performed  on  wind  tunnel  model  surface  thermocouple  data  at  the  Navy's  Hypervelocity  Wind 
Tunnel  No.  9  (Tunnel  9).  These  calculations  have  been  used  on  Tunnel  9  data  since  the  early 
1980s  when  coaxial  surface  thermocouples  began  to  replace  gardon  gauges  as  the  standard  means 
of  obtaining  heat-flux  data.  A  one-dimensional  heat  flow  assumption  allows  practical 
computation  of  heat-flux  from  a  single  thermocouple  output.  However,  this  assumption  will 
break  down  under  two-  and  three-dimensional  heat  conduction.  Quantitative  information  on  the 
limits  of  the  one-dimensional  assumption  is  needed  by  project  engineers  planning  wind  tunnel 
tests.  This  report  addresses  issues  associated  with  the  one-dimensional  assumption  and  should  be 
consulted  prior  to  making  heat  transfer  measurements  in  Tunnel  9.  Analysis  on  other  aspects  of 
the  heat  transfer  testing  in  Tunnel  9  is  ongoing. 

The  authors  acknowledge  Michael  Metzger  and  Leonard  Zentz  of  the  Weapons  Dynamics 
Division  for  their  helpful  instructions  during  the  execution  of  the  ABAQUS  finite  element  code. 
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R.  L.  SCHMIDT,  Head 

Strategic  and  Space  Systems  Department 
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INTRODUCTION 


Heat  transfer  measurements  are  a  major  component  of  the  testing  conducted  in  the  Navy's 
Hypervelocity  Wind  Tunnel  No.  9  (Tunnel  9)  located  at  the  White  Oak,  MD  site  of  the  Dahlgren 
Division,  Naval  Surface  Warfare  Center.  Measurements  made  in  Tunnel  9  have  employed  the  use 
of  coaxial  surface  thermocouples  since  the  early  1980s.*  These  transducers  are  commonly  used 
for  obtaining  transient  surface  heat-flux  measurements  on  wind  tunnel  models.^ 

Details  on  the  use  of  the  coaxial  thermocouples  are  found  in  References  1  and  2. 
Advantages  of  the  coaxial  surface  thermocouple  technique,  copied  in  part  from  Reference  2, 
include  the  following; 

a)  durability 

b)  small  size  (0.03 1 "  or  0.06 1 "  typical  diameters  available) 

c)  easy  installation,  flush  mounting  to  virtually  any  model  contour 

d)  fast  response  time  ( ^ 50  psec) 

e)  no  power  supply  required 

f)  no  calibration  required 

The  main  drawback  of  the  coaxial  thermocouple  technique  is  that  heat-flux,  the  desired  quantity, 
is  calculated  from  the  time  history  of  the  coaxial  thermocouples  output  temperature.  A  model 
must  be  constructed  in  order  to  calculate  the  heat-flux.  In  the  majority  of  Tunnel  9  tests,  good 
heat-flux  results  are  obtained  from  coaxial  thermocouples  with  relatively  simple  models. 

However,  a  single  simple  model  will  not  work  under  all  conditions.  Leading  edges,  nose  tips,  and 
other  complex  geometries  as  well  as  spatially  varying  heating  rates  can  require  different  solution 
techniques.  Understanding  the  solution  technique  and  its  assumptions  and  limitations  is  very 
important. 

The  two  common  techniques  used  to  calculate  heat -flux  from  thermocouple  outputs  are 
the  semi-infinite  slab  assumption  (see  Reference  1),  and  a  finite  thickness  one-dimensional  (ID) 
finite  difference  approximation.^’*  The  finite  difference  approximation  assumes  the  model  wall 
thickness  is  finite  and  allows  the  user  to  define  a  back-face  boundary  condition.  Since  the  finite 
difference  technique  is  more  adaptive  to  wind  tunnel  model  configurations,  it  is  the  standard 
technique  used  for  heat -flux  measurements  made  in  Tunnel  9.  QCALC  is  the  name  of  the 
subroutine  which  contains  the  finite  difference  scheme  used  for  Tunnel  9  data. 
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QCALC  SUBROUTINE 


QCALC  is  a  FORTRAN  subroutine  which  solves  the  transient  ID  heat  equation  in 
Cartesian  coordinates  (Equation  1).  Equation  2  shows  the  second  order  Euler-explicit  finite 
difference  approximation  used  in  QCALC.  Subscripts  i  and  j  refer  to  the  time  and  space  steps 
respectively.  Temperature  versus  time  data  from  the  model  surface  are  input  into  QCALC  which 
uses  Equation  2  to  solve  for  the  temperature  distribution  at  each  node  through  the  model  wall 
thickness  for  each  time  step.  Heat-flux  is  obtained  from  a  second  order  approximation  to  the 
derivative  of  the  temperature  profile  evaluated  at  the  mode!  surface. 


at  “ 
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The  solution  begins  at  time  zero  (i=l)  with  all  nodes  set  to  a  uniform  initial  temperature. 
QCALC  solves  for  the  temperature  at  successive  time  steps  by  solving  for  7;+,  j  from  Equation  2. 
The  outer  surface  boundary  condition  is  the  measured  surface  temperature.  The  back-face 
temperature  is  calculated  by  assuming  an  insulated  condition.  As  an  option,  the  back-face 
temperature  can  be  measured.  Heat-flux  is  computed  at  the  completion  of  each  time  step. 
Appendix  A  contains  the  QCALC  subroutine  FORTRAN  code.  Major  assumptions  and 
techniques  used  in  the  QCALC  model  are  summarized  below. 

a)  one-dimensional  heat  conduction  equation  in  Cartesian  coordinates 

b)  second  order  Euler-explicit  finite  difference  approximation 

c)  constant  and  homogeneous  material  properties 

d)  uniform  initial  temperature 

e)  front-face  boundary  condition  measured  by  coaxial  thermocouple 

f)  back-face  boundary  assumed  insulated  (optional  temperature  measurement) 

g)  heat-flux  is  obtained  from  a  second  order  approximation  to  the  temperature 
gradient  at  the  model  surface 

h)  negligible  radiation  effects 

Figure  1 A  illustrates  a  wind  tunnel  model  with  a  conical  geometry  and  a  spatially  varying 
heat-flux.  The  local  approximation  made  by  QCALC  is  illustrated  in  Figure  IB.  Model  surface 
curvature  and  spatial  variations  in  heat-flux  are  assumed  negligible  for  the  small  measurement 
region.  This  assumption  will  begin  to  break  down  when  the  local  spatial  variation  of  the  heating 
rate  becomes  substantial  and/or  the  geometry  of  the  model  is  not  locally  flat.  Evaluating  the 
limiting  factors  of  these  two  assumptions  is  the  objective  of  this  report. 
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Node  spacing  and  the  time  step  (^t)  must  be  chosen  based  on  the  stability  criteria*  for  this 
technique*.  The  error  in  the  second  order  technique  used  in  QCALC  is  proportional  to  the  node 
spacing  squared  (ax^)  and  can  be  improved’  by  decreasing  the  value  of  ax.  Using  too  few  nodes 
increases  ax  and  will  lead  to  larger  errors.  For  this  report  the  authors  used  a  3/8  inch  model  wall 
with  50  nodes  and  a  time  step  of  1/500  second.  The  thermal  properties  of  constantan  are  used  for 
this  analysis  and  can  be  found  in  the  QCALC  code  in  Appendix  A. 


INVESTIGATION  APPROACH 


The  approach  for  this  investigation  involved  testing  the  assumption  of  ID  conduction  in  a 
Cartesian  coordinate  system  when  applied  under  two-  and  three-dimensional  conditions.  Two 
possible  conditions  are  considered;  what  happens  when  the  surface  is  not  locally  flat,  and  what 
happens  when  the  heating  rate  varies  spatially  to  the  extent  that  it  cannot  be  considered  locally 
uniform.  The  authors  considered  various  geometries  subjected  to  spatially  uniform  heating  loads 
to  test  the  locally  flat  condition.  Spatially  varying  heating  loads  applied  to  flat  plates  were 
analyzed*  briefly  in  support  of  Tunnel  9  testing  and  those  results  are  included.  All  analysis 
assumed  homogeneous  constant  property  materials. 


MODEL  GEOMETRIES 

Typical  models  tested  in  Tunnel  9  include  many  varieties  of  cones  as  well  as  hypersonic 
airframes.  Most  of  the  surfaces  are  not  flat.  Three  geometries  were  identified  for  evaluation 
which  locally  represent  many  wind  tunnel  model  surfaces  better  than  the  flat  plate  while  remaining 
simple  enough  for  the  present  study.  The  geometries  selected  for  evaluation  were  several  sizes  of 
cones,  cylinders,  and  spheres. 


Conical  Section 


The  conical  section  is  a  small  element  of  a  seven  degree  half-angle  cone.  Directly 
applicable  to  cone  models,  this  geometry  could  also  represent  many  areas  on  the  fuselage  of 
hypersonic  aircraft.  Figure  2A  illustrates  a  conical  element.  Radius  of  curvature  and  temperature 
values  are  referenced  to  the  center  of  the  element. 


For  stability.  At  ^  (Ax)V(2  a)  where  a  is  the  thermal  diffusivity. 
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Cylindrical  Section 

The  cylindrical  section  is  simply  a  portion  of  a  circular  cylinder.  This  model  is 
representative  of  many  local  sections  of  an  airframe  fiaselage  as  well  as  rounded  leading  edge 
geometries.  Figure  2B  illustrates  a  cylindrical  element. 


Spherical  Section 

A  spherical  section  is  illustrated  in  Figure  2C.  This  geometry  is  typical  of  blunted  nose 
tips  but  could  represent  other  areas  of  wind  tunnel  models  as  well. 


PROCEDURE 

The  analysis  tested  QCALC's  ability  to  compute  heat-flux  applied  to  the  geometries 
illustrated  in  Figure  2.  Temperature  rise  data  input  to  QCALC  were  obtained  from  finite  element 
solutions*  with  known  heat-flux  inputs.  Heat-flux  calculated  by  QCALC  was  compared  to  the 
heat-flux  used  for  inputs  to  the  finite  element  solutions.  This  numerical  approach  was  chosen  to 
isolate  the  effects  of  the  ID  heat  transfer  assumptions  in  the  computation  of  heat-flux  on  complex 
models.  No  experimental  measurements  were  attempted.  The  analysis  steps  are  summarized 
below. 


Generation  of  Temperature  Data 

Temperature  versus  time  data  were  generated  using  the  ABAQUS  finite  element  code. 
For  each  case,  one  of  the  model  geometries  described  above  was  subjected  to  a  known  heat-flux 
which  varied  with  time.  The  inner  wall  (back  face)  was  insulated  (heat-flux  =  0)  to  coincide  with 
the  assumption  used  in  the  QCALC  code.  Several  finite  element  solutions  of  various  element 
sizes  and  time  steps  were  obtained  for  each  case  to  ensure  convergence  of  the  solution.  Once  a 
converged  solution  was  generated,  the  surface  temperature  at  each  time  step  was  put  into  a 
computer  file  for  subsequent  use  by  the  QCALC  code.  The  applied  heat-flux  at  each  time  step 
was  also  stored  in  a  computer  file  for  later  comparison  with  the  QCALC  results. 


Computation  of  Heat-flux 

The  temperature  histories  generated  by  the  finite  element  code  were  used  as  input  to  the 
QCALC  code.  QCALC,  as  described  above,  calculated  the  surface  heat-flux  for  each  time  step. 
The  insulated  back  wall  condition  was  used  in  QCALC  to  coincide  with  the  finite  element 


ABAQUS  finite  element  code  Version  4.9,  Hibbit,  Karlson  &  Sorenson,  Inc. 
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solutions'  back-face  boundary  condition.  Heat-flux  obtained  from  QCALC  was  stored  in  a 
computer  file  for  comparison  to  the  heat-flux  used  to  generate  the  temperature  histories  above. 


QCALC  RESULTS 


SPATIALLY  UNIFORM  HEATING  ON  VARIOUS  GEOMETRIES 

Each  of  the  geometries  illustrated  in  Figure  2  was  used  for  analyzing  the  QCALC  code. 
The  flat  plate  geometry  was  used  as  a  reference  condition.  The  heating  loads  applied  were  time 
dependent  and  spatially  uniform  over  the  modeled  geometry. 


Flat  Plate 


The  technique  was  applied  to  a  flat  plate  geometry  to  coincide  with  the  assumptions  in  the 
QCALC  code.  Results  are  shown  in  Figures  3  and  4  for  a  step  and  a  half  sine  wave  input  of  heat- 
flux.  QCALC  accurately  predicts  the  heat-flux  in  the  case  of  the  flat  plate  model. 

Figure  3  illustrates  the  normalized  output  from  QCALC  when  the  temperature  distribution 
arises  from  a  flat  plate  model  subjected  to  a  step  input  of  heat-flux.  Results  were  independent  of 
the  magnitude  of  the  step  input.  A  finite  response  time  can  be  seen  as  QCALC  responds  to  the 
step  input.  Only  a  small  deviation  in  the  calculated  heat-flux  can  be  seen  during  the  first  50 
milliseconds  of  Figure  3.  Within  10  milliseconds  from  the  initiation  of  the  step,  the  QCALC 
output  value  is  within  3  percent  of  the  expected  value  and  after  25  milliseconds  the  output  is 
within  1  percent.  These  deviations  are  considered  small  in  comparison  to  typical  heat  transfer 
uncertainties  which  are  greater  than  6  percent  (Reference  1).  It  should  be  noted  that  a  step 
increase  of  this  type  is  not  representative  of  Tunnel  9  testing.  In  addition,  the  node  spacing  and 
sampling  period  could  be  decreased  to  improve  these  results  (Reference  5). 

Figure  4  illustrates  QCALC's  ability  to  follow  a  smoothly  varying  heat-flux.  This  1/2- 
cycle-per-second  sine  wave  variation  in  heat-flux  with  a  peak  at  25  BTU/ftVsec  is  typical  of  the 
time  scale  and  magnitude  of  the  heat-flux  obtained  during  many  Tunnel  9  tests.  In  this  case, 
QCALC  predicted  the  applied  heat-flux  to  within  0.12  BTU/ftVsec. 

The  preceding  two  examples  of  QCALC  results  from  flat  plate  models  illustrate  the  ability 
of  the  QCALC  code  to  accurately  predict  heat-flux  when  the  geometry  is  locally  flat. 
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Cylindrical  Section 

Cylindrical  geometries  (Figure  2B)  of  radius  1/2, 1,  2,  4,  and  8  inches  were  subjected  to 
spatially  uniform  heating  loads.  A  step  and  a  half  sine  input  of  heat-flux  were  applied  in  the  same 
manner  as  the  flat  plate  analysis. 

Figure  5  illustrates  the  normalized  output  from  QCALC  diverging  from  the  flat  plate  result 
when  a  step  input  of  heat-flux  is  applied  to  each  of  the  investigated  model  geometries.  The  flat 
plate  results  (radius  =  «>)  are  included  from  Figure  3  as  a  reference.  The  over-prediction  of  heat- 
flux  increases  with  time  and  decreasing  radius  of  curvature. 

These  results  can  be  explained  by  comparing  the  cylindrical  geometry  to  the  flat  plate 
geometry  assumed  by  QCALC.  For  the  flat  plate,  heat  flows  perpendicular  to  the  surface  with 
constant  temperature  lines  parallel  to  the  surface.  The  cross  section  through  which  the  heat  flows 
remains  unchanged.  The  cylindrical  geometry  tends  to  funnel  the  heat  flow  into  an  ever 
decreasing  cross  sectional  area  as  the  heat  flows  inward  from  the  surface.  This  decreasing  cross 
sectional  area  makes  this  geometry  more  resistive  to  conducting  heat  away  from  the  surface  and 
hence  a  larger  surface  temperature  rise  is  observed  (when  compared  to  the  flat  plate).  The  larger 
surface  temperature  rise  will  give  rise  to  a  larger  calculated  heat-flux  from  QCALC.  For 
increased  heating  times,  the  heat  penetrates  deeper  into  the  model  and  the  problem  becomes  more 
pronounced. 

Figure  6  shows  results  for  a  half  sine  input  of  heat-flux  applied  to  several  cylindrical 
elements  of  various  geometries.  Results  from  Figure  4  (radius  =  <»)  are  included  as  a  reference. 

It  is  noted  that  the  errors  grow  with  time  and  with  decreasing  radius  of  curvature.  These  trends 
are  similar  to  those  noted  from  Figure  5  above. 


Conical  Section 


Seven  degree  half-angle  cone  elements  (Figure  2A)  of  central  radii  1/2,  1,  2,  4,  and  8 
inches  were  subjected  to  spatially  unifot  m  heating  loads.  Figure  5  shows  the  normalized  output 
from  QCALC  diverging  from  the  expected  result.  The  seven  degree  cone  models  yielded  results 
with  negligible  differences  from  the  results  generated  from  the  cylindrical  models  of  the  same 
radius. 


Spherical  Section 

Spherical  geometries  (Figure  2C)  of  radii  1,  2,  4,  8,  and  16  inches  were  subjected  to 
spatially  uniform  heating  loads.  A  step  input  of  heat-flux  was  applied  in  the  same  manner  as  the 
cone  and  cylinder  analyses.  Figure  5  shows  the  normalized  output  from  QCALC  for  the  spherical 
sections  subjected  to  a  step  input  of  heat-flux.  The  results  from  a  spherical  section  of  a  given 


6 


NSWCDD/TR-94/1 14 


radius  are  equivalent  to  results  from  a  cylindrical  section  of  half  the  radius.  For  a  given  radius, 
heat-flux  results  for  the  spherical  sections  diverge  faster  from  the  expected  results  than  the  heat- 
flux  data  from  the  cone/cylinder  sections.  The  faster  divergence  can  be  explained  by  comparing 
the  spherical  geometry  to  the  cylindrical  geometry.  As  heat  moves  inward,  away  from  the 
surface,  the  cross  section  through  which  the  heat  travels  decreases  at  a  higher  rate  for  the 
spherical  elements  when  compared  to  cylindrical  elements  of  the  same  radius.  This  will  create  a 
higher  surface  temperature  rise  for  the  spherical  element.  A  higher  surface  temperature  rise  will 
cause  a  higher  predicted  heat-flux  from  QCALC. 


SPATIALLY  VARYING  HEAT-FLUX  ON  A  FLAT  PLATE 

The  QCALC  code  was  tested  under  the  conditions  of  a  spatially  varying  heat-flux 
(Reference  6).  Two-dimensional  (2D)  solutions  were  generated  for  a  heat-flux  which  varied 
linearly  in  space  along  the  face  of  a  flat  plate  element.  The  QCALC  code  was  used  to  compute 
the  heating  load  at  one  point  on  the  surface.  Deviations  of  the  ID  QCALC  solution  compared  to 
the  2D  solution  are  given  in  Figure  7  which  was  obtained  from  Reference  6.  The  error  is  shown 
to  increase  with  time  and  with  the  magnitude  of  the  spatial  gradient  of  heat-flux.  These  results 
can  be  used  to  show  the  reliability  of  the  ID  heat-flux  calculation  when  used  for  spatially  varying 
heating  loads.  As  shown  in  Figure  7,  extreme  spatial  variations  in  the  heat-flux  will  lead  to  errors 
in  the  calculated  heat-flux.  Typical  values  of  (dQ/ds)  /  Q*  for  Tunnel  9  testing  are  much  less  than 
1. 


QCALC  MODIFICATIONS 


The  QCALC  code  works  well  when  used  on  geometries  which  locally  can  be  considered 
flat  and  when  the  spatial  gradient  of  the  heat-flux  is  not  too  large.  It  has  been  shown  in  the 
previous  section  that  the  QCALC  code  will  diverge  from  the  expected  results  when  these 
conditions  are  not  met.  Modified  codes  have  been  developed  to  solve  the  ID  heat  equation  in 
cylindrical  and  spherical  coordinates  under  spatially  uniform  heating  loads.  These  codes  are 
useful  when  model  geometry  can  locally  be  represented  by  cylindrical/conical  or  spherical 
elements. 


*  Q  is  the  local  heat-flux  value  and  s  is  the  direction  along  the  model  surface. 
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CYLINDRICAL  COORDINATES  (QCYL) 

For  cases  where  model  geometry  is  best  described  as  a  cone  or  cylinder,  the  QCALC  code 
was  modified  to  cylindrical  coordinates.  The  governing  ID  Cartesian  heat  equation  (Equation  1) 
was  replaced  with  its  ID  representation  in  cylindrical  coordinates  (Equation  3).  A  second  order 
Euler-explicit  finite  difference  scheme  was  used  to  approximate  Equation  3 .  This  new  code  is 
referred  to  as  QCYL  and  is  available  for  use  in  support  of  Tunnel  9  testing.  Modifications  made 
to  the  QCALC  code  to  create  the  QCYL  code  are  contained  in  Appendix  B. 


^  -  L  iL 

r  dr  U  dt 


(3) 


The  QCYL  program  is  based  on  cylindrical  model  surfaces  which  locally  represent  many 
of  the  models  tested  in  Tunnel  9.  A  comparison  of  QCALC  results  with  QCYL  results  is  shown 
in  Figure  8  for  a  cylindrical  shell  model  with  1"  radius  subjected  to  a  step  input  of  heat -flux. 
Results  for  several  radii  were  tested  with  similar  results.  The  QCYL  program  takes  into  account 
the  cylindrical  geometry  of  the  element  and  does  not  over-predict  the  heat-flux.  As  mentioned 
earlier,  the  conical  and  cylindrical  geometries  yield  results  which  are  essentially  the  same.  The 
QCYL  program  is  therefore  applicable  to  cases  where  model  geometry  can  be  described  as  a 
cone*  or  cylinder. 

The  response  of  the  QCYL  program  to  a  step  input  of  heat-flux  on  a  cylindrical  model  is 
slightly  slower  than  the  response  of  QCALC  to  a  step  input  of  heat-flux  on  a  flat  plate.  Figure  9 
shows  the  results  of  QCALC  and  QCYL  for  a  step  input  of  heat-flux  on  a  flat  plate  and  1  inch 
radius  cylindrical  shell  respectively.  Both  calculations  were  executed  with  the  same  time  step  and 
node  spacing.  The  QCALC  code  responds  slightly  faster  to  the  step  input  when  compared  to  the 
cylindrical  QCYL  program.  It  is  assumed  by  the  authors  that  the  extra  term  in  Equation  3  (1/r 
dT/dr)  adds  some  numerical  dissipation  to  the  QCYL  solution.  This  factor  should  be  considered 
when  using  the  QCYL  program  to  measure  rapid  changes  in  heating  loads.  The  response  shown 
in  Figure  9  is  considered  adequate  for  most  Tunnel  9  testing  since  typical  model  surface  heat-flux 
values  vary  much  slower  than  the  evaluated  step  input.  Increasing  the  number  of  nodes  used  as 
well  as  the  sampling  frequency  would  allow  the  QCYL  program  to  respond  faster. 


SPHERICAL  COORDINATES  (QSPH) 

Next  the  QCALC  code  was  modified  to  allow  the  input  of  radius  of  curvature  for  the 
spherical  element  case.  The  governing  ID  Cartesian  heat  equation  (Equation  1)  was  replaced 
with  its  representation  in  spherical  coordinates  (Equation  4).  A  second  order  Euler-explicit  finite 


Data  from  a  seven  degree  half-angle  cone.  Larger  angles  would  require  further  study. 
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difference  scheme  was  used  to  approximate  Equation  4.  This  new  code  is  referred  to  as  QSPH 
and  is  available  for  use  in  support  of  Tunnel  9  testing.  Modifications  made  to  the  QCALC  code  in 
the  development  of  the  QSPH  code  are  contained  in  Appendix  C. 


^  .  L  iL 

dr^  r  Br  a  dt 

It  was  noted  from  Figure  5  that  a  sphere  radius  twice  as  large  as  a  given  cylindrical  radius 
will  have  the  same  effect  on  the  ID  QCALC  results.  Equation  4  shows  the  origins  of  this 
relationship.  The  second  term  in  Equation  4  is  simply  twice  the  second  term  xfom  the  cylindrical 
coordinates  equation  (Equation  3).  This  second  term  is  twice  as  important  in  the  equation  for 
spherical  geometries. 

Figure  10  illustrates  the  ability  of  the  QSPH  code  to  predict  the  heat-flux  applied  to  a 
spherical  element  with  1  inch  radius.  The  data  in  Figure  10  are  not  as  smooth  looking  as  the  data 
presented  earlier.  This  result  is  due  to  the  numerical  format  of  the  input  temperature  data  for  this 
individual  case  and  not  the  result  of  the  QSPH  program.  Original  input  data  were  deleted  and 
only  a  copy  of  the  input  temperature  data,  truncated  to  three  decimal  positions,  was  available  for 
analysis.  This  truncation  leads  to  the  slight  variations  in  the  heat-flux  since  it  is  obtained  through 
taking  derivatives. 


CONCLUSIONS 


The  validity  of  the  assumptions  in  the  QCALC  code  should  be  considered  when  analyzing 
any  heat-flux  data  from  Tunnel  9.  In  particular,  the  ID  conduction  (Cartesian)  assumption  begins 
to  break  down  when  the  model  surface  curvature  becomes  too  small  and/or  the  heating  load 
spatial  gradient  becomes  too  large.  Longer  run  times  lead  to  greater  deviations  in  both  of  these 
cases.  The  QCALC  program  should  be  used  for  most  of  the  measurements  made  in  Tunnel  9. 
However,  the  QCYL  (or  QSPH)  program  will  do  a  better  job  of  computing  heat-flux  on  model 
configurations  when  the  local  cylindrical  (or  spherical)  radius  of  curvature  is  small.  Spherical 
effects  become  important  at  radii  twice  as  large  as  the  radii  where  cylindrical  effects  become 
important.  It  is  also  noted  that  more  nodes  would  be  required  in  the  QSPH  and  QCYL  programs 
to  obtain  the  response  time  characteristics  of  the  QCALC  code. 
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RECOMMENDATIONS 


When  considering  results  in  this  report,  the  method  of  analysis  must  be  considered.  One 
idealized  model  is  being  compared  to  a  second  idealized  model.  No  physical  measurements  were 
made. 


The  finite  element  solutions  modeled  homogeneous  material  sections.  Most  of  the  surface 
temperature  measurements  made  in  Tunnel  9  involve  the  use  of  coaxial  thermocouples.  For  this 
type  of  gauge,  the  measurement  section  is  not  homogeneous.  A  typical  coaxial  gauge  consists  of 
a  0.012-inch  diameter  constantan  wire  with  a  0.0005-inch  thick  insulative  coating  protecting  it 
from  the  0.0607-inch  diameter  chrome!  outer  jacket.  This  gauge  is  cemented  into  a  0.0625  inch 
hole  in  stainless  steel.  The  effects  of  having  three  different  metals,  the  cement,  and  possible 
contact  resistances  between  the  metals  were  not  considered  in  this  report.  Further  study  of  these 
factors  along  with  experimental  measurements  are  recommended  to  further  understand  the  heat 
transfer  results  obtained  from  coaxial  thermocouple  gauges. 
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(A)  TYPICAL  WIND  TUNNEL  MODEL  WITH  VARIABLE  HEATING 


q(t)  heating  rate  varies  with  time 


local  s| 

patially  uniform  heating  1 

load 
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Conditions  for  ID  heat  equation  solution 
Initial  Condition: 

T(x,t)  =  initial  temperature  (constant  and  uniform) 

Boundary  Conditions: 

X  =  0 

T  =  T  measured,  surface  thermocouple 

X  =  d 

dT/dx  =  0,  assumed  insulated 
or 

T  =  T  initial,  assumed  constant 
or 

T  =  T  measured,  backface  thermocouple 


thermocouple  output 


(B)  IDEALIZED  LOCAL  GEOMETRY  AND  BOUNDARY  CONDITIONS 

FIGURE  1.  ILLUSTRATION  OF  LOCAL  APPROXIMATION  MADE  BY  QCALC 
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conical  element 


(A)  CONICAL  SECTION 


(B)  CYLINDRICAL  SECTION 


(C)  SPHERICAL  SECTION 


FIGURE  2.  GEOMETRIC  SECTIONS  USED  FOR  ANALYSIS 


Q  (BTUmscc) 


FIGURE  3.  NORMALIZED  QCALC  OUTPUT  FOR  STEP  INPUT  APPLIED  TO  FLAT  PLATE 


FIGURE  4.  QCALC  COMPARISON  FOR  SINUSOIDAL  HEAT-FLUX  APPLIED  TO  FLAT  PLATE 


Qcalc/Q, 
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nCURE  5.  NORMALIZED  QCALC  OUTPUT  FOR  STEP  INPUT  APPLIED  TO  MODEL  GEOMETRIES 
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FIGURE  6.  QCALC  RESULTS  FOR  SINUSOIDAL  HEAT-FLUX  APPLIED  TO  CYLINDRICAL  MODELS 
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FIGURE  8,  COMPARISON  OF  NORMALIZED  QCALC  AND  QCYL  RESULTS 
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APPENDIX  A 

QCALC  FORTRAN  CODE 


Q********************* ********************************************** ******* 

SUBROUTINE  QCALC (TIME,  Tl,  T2,  QOUT,  NTHI,  NODES, THICK) 

CCCCC  call  qcalc  (tim,  CO  (1,  k)  ,  bface,  qout,  nptq,  nodes,  thiclcft) 

C 

C  Now  the  standard  version...  2/88 

C 

C  PURPOSE 

C  ======= 

C  THIS  SUBROUTINE  USES  THE  FINITE  (FORWARD)  DIFFERENCE  METHOD  TO 

C  DETERMINE  SURFACE  HEAT  TRANSFER  FROM  THE  CALCULATED  TEMPERATURE  PROFILE. 

C  IT  STARTS  AT  SECOND  POINT  IN  THE  T  ARRAYS  AND  CALCULATES  THE  HEAT 

C  FLUX, QDOT,  FOR  NTHI  PTS . 

C 

C  EXPLANATION 

C  =========== 

C  STABILITY  CONSTRAINTS  OF  THE  FINITE  DIFFERENCE  METHOD  MUST  BE  SATISFIED 

C  (THET  <=  .5,  ie.  THE  TIME  INTERVAL  MUST  BE  SMALL  ENOUGH), 

C  IF  THE  RECORDED  TIME  INTERVAL  IS  TOO  LARGE,  THEN  EACH 

C  INTERVAL  IS  DIVIDED  INTO  SUBTINTERVALS . 

C 

C  FORWARD  FINITE  DIFFERENCE  DIAGRAM 

C 

C  t :  time 

C  T:  temperature 

C  x:  distance 

C 
C 

C  t(j)  !  T4 

C  t(j-l)  !  Tl  T2  T3 

C  ! 

C  - 

C  x(i-l)  x(i)  x(i+l) 

C 

C  temperatures  T1,T2,T3  are  used  to  calculate  the 

C  temperature  T4 

C 

C  SUBROUTINES  CALLED 

C  (NONE) 

C 

C  VARIABLE  GLOSSARY 

C  ================= 

C  ADJUST  =  AMOUNT  TO  SUBTRACT  FROM  THE  BACKSIDE  THERMOCOUPLE 

C  TEMPERATURES  IN  'T2' 

C  COND  =  CODUCTIVITY  (BTU/FT-SEC-DEGF) 

C  CONDODX  =  CONDUCTIVITY/DX 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


- PRELIMINARY  CALCULATIONS 

IF(NODES.EQ.O)  NODES  =  10 
IF (NODES, LE.MAXNODE)  GO  TO  10 

PRINT  *,  '  SUBROUTINE  QCALC...  NUMBER  OF  NODES  GREATER  THAN  DIMEN 

ISION;  NUMBER  OF  NODES  SET  EQUAL  TO',  MAXNODE 
NODES  =  MAXNODE 
10  NODEMl  =  NODES  -  1 

HIDT  =  TIME (2)  -  TIME(l) 

DX  =  THICK  /  NODEMl 

assume  properties  of  steel  model,  and  chromel-constantan  TCP's  are  close 
enough  to  use  constantan  properties 

These  are  the  CONSTANTAN  properties  used  for  data  reduction 
CP=.094  ! (BTU/Lbra  F) 

RHO=.322*172r  !  (Lbm/in^3)  *  1728  in''3/ft''3 

COND=(0.2676E-03) *12  ! (BTU/in  sec  F)  *  12  in/ft 

CONODX  =  COND  /  DX 

THETOT  =  COND  /  (CP  *  RHO  *  DX**2) 

-  THET  MUST  BE  LESS  THAN  0.50  TO  MAKE  SOLUTION  PROCEDURE  STABLE 

THET  =0.5 

DTLO  =  THET  /  THETOT 
NTLO  =  HIDT  /  DTLO 
IF  (NTLO  .LT,  0)  THEN 


CP 

DTLO 

DX 

HIDT 

IHI 

INODE 

ITLO 

LI 

L2 

LHOLD 

NODEMl 

NODES 

NTLO 

QOUT 

RHO 

T 


T1 

T2 

TFAC 

THET 

THETOT 

THICK 

TIME 


SPECIFIC  HEAT  (BTU/LB-DEGF) 

DELTA  TIME  FOR  A  SUBINTERVAL 

LENGTH  BETWEEN  TWO  NODES 

NUMBER  OF  RECORDED  TIMES 

COUNTER  OF  TIME  INTERVALS 

COUNTER  OF  NODES 

COUNTER  OF  TIM..  SUBINTERVALS 

INDEX  OF  TEMPERATURE  T,  REPRESENTS  t(j-l) 

INDEX  OF  TEMPERATURE  T,  REPRESENTS  t(j) 

TEMPORARY  STORAGE  USED  TO  INTERCHAGE  LI  AND  L2 
NODES  -  1 

NUMBER  OF  X  POINTS 

NUMBER  OF  SUBINTERVALS  IN  EACH  TIME  INTERVAL 
HEAT  FLUX  (BTU/FT2-SEC) 

DENSITY  (LB/FT3) 

ARRAY  WITH  EACH  ELEMENT  CONTAINING  THE  TEMPERATURE 
AT  A  SPECIFIC  LOCATION  AND  TIME.  ONLY  THE 
TEMPERATURES  AT  TIMES  t(j)  and  t(j-l)  ARE  STORED 
ARRAY  OF  COAX  GAGE  TEMPERATURES  AT  THE  TIMES  CONTAINED 
IN  'TIME' 

ARRAY  OF  BACKSIDE  THERMOCOUPLE  TEMPERATURES  AT  THE 
TIMES  CONTAINED  IN  'TIME' 

FRACTION  OF  TIME  INTERVAL  USED  TO  PRODUCE  THE 
TEMPERATURE  AT  A  SUBINTERVAL 
(K  *  dT)  /  (RHO  *  C  *  dx) 

K  /  (RHO  *  C  *  dx**2) 

LENGTH  OF  SLAB  (ft) 

ARRAY  OF  RECORDED  INPUT  TIMES 


PARAMETER (  MAXNODE=50  ) 

DIMENSION  TIME(l) ,T1 (1) ,T2 (1) ,QOUT(l) , T (MAXNODE,  2 ) 
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NTLO  =  lABS(NTLO) 

ENDIF 

NTLO  =  NTLO  +  1 
DTLO  =  HIDT  /  NTLO 
THET  =  THETOT  *  DTLO 
B  =  1.0  -  (2.0  *  THET) 

C 

C  if  semi-infinite  slab  is  assumed,  then  in  main  program  set  T2  array  of 

C  backface  temps  equal  to  first  point  of  T1  array  of  frontface  temps 

ADJUST  =  T2(l)  -  Tl(l)  !  backface (1)  -  co(l,k) 

DO  52  I  =  1,NTHI  !  I=l,nptq 

T2(I)  =  T2(I)  -  ADJUST  !  backface  temps 
52  CONTINUE 
C 

DO  102  I  =  1, NODES 

T(I,1)  =  Tl(l)  !  co(l,k)  thru  all  nodes 

102  CONTINUE 
C 

LI  =  1 
L2  =  2 
C 

C - HI  TIME  LOOP  (loop  through  time  intervals) 

C 

DO  302  IHI  =  2,NTHI  !  Ihi=l,nptq 

C 

C - LO  TIME  LOOP  (loop  through  time  subtintervals ) 

C 

DO  252  ITLO  =  l,NTLO 
C 

C - CALCULATE  TEMPERATURE  FOR  CURRENT  TIME  AT  ALL  NODES 

DO  202  INODE  =  2,N0DEM1 

T (INODE, L2)  =  THET  *  (T (INODE  +  1 , LI )  +  T (INODE-1,  LI )  ) 

+  +  B  *  T (INODE, LI) 

202  CONTINUE 

TFAC  =  FLOAT (ITLO)  /  FLOAT (NTLO) 

T(1,L2)  =  Tl(IHI-l)  +  (Tl(IHI)  -  Tl(IHI-l))  *  TFAC 
C 

C  Longer  times  sc  can't  assume  semi-infinite  slab  equation 

C  Qdot  =  0  at  backface  rather  than  a  constant 

T(NODES,L2)  =  T (NODES-1 , L2) - (T (NODES-2 , L2 ) -T (NODES-1, L2 ) ) / 3 .  !  SD5 

CCCCCCCCCCC  T (NODES, L2)=  T2(IHI-1)  +  (T2(IHI)  -  T2(IHI-1))  *  TFAC  !  semi-inf 


C 

LHOLD  =  LI 
LI  =  L2 
L2  =  LHOLD 
252  CONTINUE 

CCCCCC 

C - NEW  WALL  TEMPERATURES  NOW  AT  LI  (ie.  at  the  input  recorded  time, 

C - TIME  (IHI)  ) 

C 

CCCCCCCC  QOUT(IHI)  =  CONODX  *  (T(1,L1)  -  T(2,L1)) 

CCCCCCCC  QOUT(IHI)  =  CONODX  *  (-2.  *  T(1,L1)  +  3.  *  T(2,L1) 

C  +  -  T(3,L1)  ) 

C 

QOUT  (IHI ) =-l . * (CONODX/2 . ) * (-3 . *T { 1, LI ) +4 . *T (2, LI) 

+  -T(3,L1)) 


302  CONTINUE 
RETURN 
END 
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APPENDIX  B 

QCYL  FORTRAN  CODE  MODIFICATIONS 


SUBROUTINE  QCYL (TIME,  Tl,  T2,  QOUT,  NTHI,  NODES,  THICK,  R) 


C 

C 

C 

C 

C 

C 


R 


Radius  of  curvature  of  surface  (Ft) 

Gage  is  normally  mounted  on  outside  of  model  and  local 
radii  will  decrease  from  1  to  nodes. 

If  gage  is  mounted  on  inside  of  model,  then  flag  code 
by  inputting  Radius  as  a  negative  value.  Local  radii 
will  increase  from  1  to  nodes. 


DIMENSION  TIME(l),  Tl(l),  T2(l),  QOUT(l),  T (MAXNODE, 2) , 
+  RC (MAXNODE) 


-  Calculate  the  local  radial  value  for  each  node 

If  R  is  input  as  negative,  then  curvature  is  outward 


DO  I  =  1, NODES 

DELR  =  (I-l)  *  DX 
RC(I)  =  R  -  DELR 
ENDDO 

IF(R  .LT.  0.)  THEN 
DO  I  =  1, NODES 

RC(I)  =  -  RC(I) 
ENDDO 
ENDIF 


! 

!  radius  of  curvature  vector  (ie  1,.9,.8,..., 
!  outward  curvature  ,  convex 
!  (ie  1.0,  1.1,  1.2, . . . ) 


C -  Calculate  temperature  for  current  time  at  all  nodes 

TFAC  =  FLOAT (ITLO)  /  FLOAT (NTLO) 

T(1,L2)  =  Tl(IHI-l)  +  (Tl(IHl)  -  Tl(lHI-l))  *  TFAC 
DO  INODE  =  2,NODEMl 

TPART  =  THET  *  (T (INODE+1, LI)  +  T (INODE-1,  LI) )  + 

+  B  *  T (INODE, LI) 

T(INODE,L2)  =  TPART  -  (T (INODE+l, LI )  -  T (INODE-1,  LI) )  * 

+  THET  *  DX  /  RC (INODE)  /  2.0  !  cyl  coords  ID 

ENDDO 
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APPENDIX  C 

QSPH  FORTRAN  CODE  MODIFICATIONS 


SUBROUTINE  QSPH (TIME,  Tl,  T2,  QOUT,  NTHI,  NODES,  THICK,  R) 


C 

C 

C 

C 

C 

C 


R 


Radius  of  curvature  of  surface  (Ft) 

Gage  is  normally  mounted  on  outside  of  model  and  local 
radii  will  decrease  from  1  to  nodes. 

If  gage  is  mounted  on  inside  of  model,  then  flag  code 
by  inputting  Radius  as  a  negative  value.  Local  radii 
will  increase  from  1  to  nodes. 


DIMENSION  TIME (1) ,  Tl(l),  T2(l),  QOUT(l),  T (MAXNODE, 2 ) , 
+  RC (MAXNODE) 


-  Calculate  the  local  radial  value  for  each  node 

If  R  is  input  as  negative,  then  curvature  is  outward 


DO  I  =  1, NODES 

DELR  =  (I-l)  *  DX 
RC(I)  =  R  -  DELR 
ENDDO 

IF(R  .LT.  0.)  THEN 
DO  I  =  1, NODES 

RC(I)  =  -  RC(I) 
ENDDO 
ENDIF 


!  radius  of  curvature  vector  (ie  1,.9,.8,..., 
!  outward  curvature  ,  convex 
!  (ie  1.0,  1.1,  1.2, . . .) 


C -  Calculate  temperature  for  current  time  at  all  nodes 

TFAC  =  FLOAT (ITLO)  /  FLOAT (NTLO) 

T(1,L2)  =  Tl(IHI-l)  +  (Tl(IHI)  -  Tl(IHI-l))  *  TFAC 
DO  INODE  =  2,N0DEM1 

TPART  =  THET  *  (T (INODE+1, LI)  +  T(IN0DE-1, LI) )  + 

+  B  *  T (INODE, LI) 

T(INODE,L2)  =  TPART  -  (T (INODE+1, LI )  -  T (INODE-1,  LI) )  * 
+  .  THET  *  DX  /  RC{ INODE)  !  sph  coords  ID 

ENDDO 
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